Warning: package 'terra' was built under R version 4.3.3
library(ncdf4)library(chron)library(viridis)library(readxl)library(ggsidekick)theme_set(theme_sleek()) # devtools::install_github("seananderson/ggsidekick") # not on CRAN# Source code for map plotsdevtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
Warning: package 'sf' was built under R version 4.3.3
plot_map +geom_point(data =filter(pred_grid, year ==1999), aes(X *1000, Y *1000), size =0.001, alpha =0.5) +NULL
Warning: Removed 2665 rows containing missing values or values outside the scale range
(`geom_point()`).
# Add lat and lon# Need to go from UTM to lat long for this one...# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-rxy <-as.matrix(pred_grid |> dplyr::select(X, Y) |>mutate(X = X *1000, Y = Y *1000))v <-vect(xy, crs ="+proj=utm +zone=33 +datum=WGS84 +units=m")y <-project(v, "+proj=longlat +datum=WGS84")lonlat <-geom(y)[, c("x", "y")]pred_grid$lon <- lonlat[, 1]pred_grid$lat <- lonlat[, 2]ggplot(filter(pred_grid, year ==1999), aes(lon, lat)) +geom_point()
# Add depth now to remove islands and remaining land# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package# https://emodnet.ec.europa.eu/geoviewer/dep_raster <- terra::rast(paste0(home, "/data/covariates/depth/Mean depth natural colour (with land).nc"))class(dep_raster)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
crs(dep_raster, proj =TRUE)
[1] "+proj=longlat +datum=WGS84 +no_defs"
plot(dep_raster)
pred_grid$depth <- terra::extract(dep_raster, pred_grid |> dplyr::select(lon, lat))$elevationggplot(pred_grid, aes(lon, lat, color = depth *-1)) +geom_point()
pts <-SpatialPoints(cbind(dat$lon, dat$lat),proj4string =CRS(proj4string(shape)))dat$subdiv <-over(pts, shape)$Area_27# Rename subdivisions to the more common names and do some more filtering (by sub div and area)sort(unique(dat$subdiv))
# Read data on rectangle levelspr <-read_xlsx(paste0(home, "/data/covariates/pelagic/N and B per Rect. 1991-2023.xlsx"),sheet =4) |>rename(ices_rect = RECT,year = ANNUS ) |>mutate_at(vars(`1`, `2`, `3`, `4`, `5`, `6`, `7`, `8`), ~replace_na(., 0)) |># I need to replace NA with 0, else I can't sum! According to Olavi who sent the data, NA means 0 and nothing else. Rectangle*year combinations that do not have information about biomass are simply not included in this datamutate(ices_rect =as.factor(ices_rect),Species ="Sprat",biomass_spr =`1`+`2`+`3`+`4`+`5`+`6`+`7`+`8`,id_r =paste(ices_rect, year, sep ="_") ) |>filter(year >=1993) |> dplyr::select(year, ices_rect, biomass_spr) |>summarise(biomass_spr =mean(biomass_spr), .by =c(ices_rect, year))her <-read_xlsx(paste0(home, "/data/covariates/pelagic/N and B per Rect. 1991-2023.xlsx"),sheet =3) |>rename(ices_rect = RECT,year = ANNUS ) |>mutate_at(vars(`1`, `2`, `3`, `4`, `5`, `6`, `7`, `8`), ~replace_na(., 0)) |># I need to replace NA with 0, else I can't sum! According to Olavi who sent the data, NA means 0 and nothing else. Rectangle*year combinations that do not have information about biomass are simply not included in this datamutate(ices_rect =as.factor(ices_rect),Species ="Sprat",biomass_her =`1`+`2`+`3`+`4`+`5`+`6`+`7`+`8`,id_r =paste(ices_rect, year, sep ="_") ) |>filter(year >=1993) |> dplyr::select(year, ices_rect, biomass_her) |>summarise(biomass_her =mean(biomass_her), .by =c(ices_rect, year))pelagics <-left_join(spr, her, by =c("year", "ices_rect"))dat <- dat |> tidylog::left_join(pelagics, by =c("year", "ices_rect"))# If there are no data in a specific rectangle, use the sub-division mean. If no values in the sub div, use annual meandat <- dat |>mutate(biomass_spr_sd_mean =mean(biomass_spr, na.rm =TRUE),biomass_her_sd_mean =mean(biomass_her, na.rm =TRUE),.by =c(sub_div, year) ) |>mutate(biomass_spr =ifelse(is.na(biomass_spr), biomass_spr_sd_mean, biomass_spr),biomass_her =ifelse(is.na(biomass_her), biomass_her_sd_mean, biomass_her) ) |>mutate(biomass_spr_yr_mean =mean(biomass_spr, na.rm =TRUE),biomass_her_yr_mean =mean(biomass_her, na.rm =TRUE),.by =c(year) ) |>mutate(biomass_spr =ifelse(is.na(biomass_spr), biomass_spr_yr_mean, biomass_spr),biomass_her =ifelse(is.na(biomass_her), biomass_her_yr_mean, biomass_her) )
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/covariates/oxygen/cmems_mod_bal_bgc_my_P1M-m_1718018848983.nc (NC_FORMAT_NETCDF4):
1 variables (excluding dimension variables):
float o2b[longitude,latitude,time] (Contiguous storage)
units: mmol m-3
_FillValue: NaN
standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
long_name: Dissolved Oxygen Concentration
3 dimensions:
latitude Size:404
standard_name: latitude
long_name: Latitude
units: degrees_north
unit_long: Degrees North
axis: Y
valid_min: 53.0082969665527
valid_max: 65.8909912109375
longitude Size:505
standard_name: longitude
long_name: Longitude
units: degrees_east
unit_long: Degrees East
axis: X
time Size:348
standard_name: time
long_name: Time
units: seconds since 1970-01-01 00:00:00
calendar: gregorian
axis: T
12 global attributes:
Conventions: CF-1.11
title: CMEMS ERGOM monthly integrated model fields
institution: Baltic MFC, PU Danish Meteorological Institute
source: CMEMS BAL MFC NEMO model output converted to NetCDF
history: Fri Dec 16 10:51:30 2022: cdo -z zip_1 --sortname copy /net/isilon/ifs/arch/home/iri/BALMFC/Nemo4.0/BALMFCMYP2022_univariateSST//BAL-ERGOM_BGC-DailyMeans-19930101.nc tmp_19930101.nc
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
comment: Data on cropped native product grid. Horizontal velocities destagged
subset:source: ARCO data downloaded from the Marine Data Store using the MyOcean Data Portal
subset:productId: BALTICSEA_MULTIYEAR_BGC_003_012
subset:datasetId: cmems_mod_bal_bgc_my_P1M-m_202303
subset:date: 2024-06-10T11:27:28.983Z
oxy_tibble <-tidync(paste(covPath, "oxygen","cmems_mod_bal_bgc_my_P1M-m_1718018848983.nc",sep ="/")) |>hyper_tibble() |>mutate(date =as_datetime(time, origin ="1970-01-01")) |>mutate(month =month(date),day =day(date),year =year(date),month_year =paste(month, year, sep ="_"),oxy = o2b *0.0223909 )# Now do recent data (forecast)print(nc_open(paste(covPath, "oxygen", "cmems_mod_bal_bgc_anfc_P1M-m_1718022118172.nc", sep ="/")))
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/covariates/oxygen/cmems_mod_bal_bgc_anfc_P1M-m_1718022118172.nc (NC_FORMAT_NETCDF4):
1 variables (excluding dimension variables):
float o2b[longitude,latitude,time] (Contiguous storage)
units: mmol m-3
_FillValue: NaN
standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
long_name: Dissolved Oxygen Concentration
3 dimensions:
latitude Size:323
standard_name: latitude
long_name: Latitude
units: degrees_north
unit_long: Degrees North
axis: Y
valid_min: 53.0082969665527
valid_max: 65.8909912109375
longitude Size:456
standard_name: longitude
long_name: Longitude
units: degrees_east
unit_long: Degrees East
axis: X
time Size:30
standard_name: time
long_name: Time
units: seconds since 1970-01-01 00:00:00
calendar: gregorian
axis: T
11 global attributes:
Conventions: CF-1.11
title: CMEMS ERGOM monthly integrated model fields
institution: Baltic MFC, PU Swedish Meteorological and Hydrological Institute
source: CMEMS BAL MFC NEMO model output converted to NetCDF
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
comment: Data on cropped native product grid. Horizontal velocities destagged
subset:source: ARCO data downloaded from the Marine Data Store using the MyOcean Data Portal
subset:productId: BALTICSEA_ANALYSISFORECAST_BGC_003_007
subset:datasetId: cmems_mod_bal_bgc_anfc_P1M-m_202311
subset:date: 2024-06-10T12:21:58.172Z
oxy_tibble_new <-tidync(paste(covPath, "oxygen","cmems_mod_bal_bgc_anfc_P1M-m_1718022118172.nc",sep ="/")) |>hyper_tibble() |>mutate(date =as_datetime(time, origin ="1970-01-01")) |>mutate(month =month(date),day =day(date),year =year(date),month_year =paste(month, year, sep ="_"),oxy = o2b *0.0223909 ) |>filter(year >2021)oxy_tibble <-bind_rows(oxy_tibble, oxy_tibble_new)# Loop through all year combos, extract the temperatures at the data locationsoxy_list <-list()for (i inunique(dat$month_year)) { d_sub <-filter(dat, month_year == i) oxy_tibble_sub <-filter(oxy_tibble, month_year == i)# Convert to rasterggplot(oxy_tibble_sub, aes(longitude, latitude)) +geom_point(size =0.1) oxy_raster <-as_spatraster(oxy_tibble_sub,xycols =2:3,crs ="WGS84", digits =3 )ggplot() +geom_spatraster(data = oxy_raster$oxy, aes(fill = oxy)) +scale_fill_viridis(option ="magma") +ggtitle(i)# Extract from raster d_sub$oxy <- terra::extract( oxy_raster$oxy, d_sub |> dplyr::select(lon, lat) )$oxy# Save oxy_list[[i]] <- d_sub}d_oxy <-bind_rows(oxy_list)
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/covariates/temperature_salinity/cmems_mod_bal_phy_my_P1M-m_1718087629163.nc (NC_FORMAT_NETCDF4):
2 variables (excluding dimension variables):
float bottomT[longitude,latitude,time] (Contiguous storage)
units: degrees_C
_FillValue: NaN
standard_name: sea_water_potential_temperature_at_sea_floor
long_name: Sea floor potential temperature
float sob[longitude,latitude,time] (Contiguous storage)
units: 0.001
_FillValue: NaN
standard_name: sea_water_salinity_at_sea_floor
long_name: Sea water salinity at sea floor
3 dimensions:
latitude Size:323
standard_name: latitude
long_name: Latitude
units: degrees_north
unit_long: Degrees North
axis: Y
valid_min: 53.0082969665527
valid_max: 65.8909912109375
longitude Size:418
standard_name: longitude
long_name: Longitude
units: degrees_east
unit_long: Degrees East
axis: X
time Size:348
standard_name: time
long_name: Time
units: seconds since 1970-01-01 00:00:00
calendar: gregorian
axis: T
11 global attributes:
Conventions: CF-1.11
title: CMEMS NEMO monthly integrated model fields
institution: Baltic MFC, PU Danish Meteorological Institute
source: CMEMS BAL MFC NEMO model output converted to NetCDF
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
comment: Data on cropped native product grid. Horizontal velocities destagged
subset:source: ARCO data downloaded from the Marine Data Store using the MyOcean Data Portal
subset:productId: BALTICSEA_MULTIYEAR_PHY_003_011
subset:datasetId: cmems_mod_bal_phy_my_P1M-m_202303
subset:date: 2024-06-11T06:33:49.163Z
st_tibble <-tidync(paste(covPath, "temperature_salinity","cmems_mod_bal_phy_my_P1M-m_1718087629163.nc",sep ="/")) |>hyper_tibble() |>mutate(date =as_datetime(time, origin ="1970-01-01")) |>mutate(month =month(date),day =day(date),year =year(date),month_year =paste(month, year, sep ="_") )# Now do recent data (forecast)print(nc_open(paste(covPath, "temperature_salinity", "cmems_mod_bal_phy_anfc_P1M-m_1718087127439.nc", sep ="/")))
File /Users/maxlindmark/Dropbox/Max work/R/pred-prey-overlap/data/covariates/temperature_salinity/cmems_mod_bal_phy_anfc_P1M-m_1718087127439.nc (NC_FORMAT_NETCDF4):
2 variables (excluding dimension variables):
float bottomT[longitude,latitude,time] (Contiguous storage)
units: degrees_C
_FillValue: NaN
standard_name: sea_water_potential_temperature_at_sea_floor
long_name: Sea water potential temperature at sea floor (given for depth comprise between 0 and 500m)
float sob[longitude,latitude,time] (Contiguous storage)
units: 0.001
_FillValue: NaN
standard_name: sea_water_salinity_at_sea_floor
long_name: Sea water salinity at sea floor
3 dimensions:
latitude Size:395
standard_name: latitude
long_name: Latitude
units: degrees_north
unit_long: Degrees North
axis: Y
valid_min: 53.0082969665527
valid_max: 65.8909912109375
longitude Size:426
standard_name: longitude
long_name: Longitude
units: degrees_east
unit_long: Degrees East
axis: X
time Size:30
standard_name: time
long_name: Time
units: seconds since 1970-01-01 00:00:00
calendar: gregorian
axis: T
11 global attributes:
Conventions: CF-1.11
title: CMEMS NEMO monthly integrated model fields
institution: Baltic MFC, PU Swedish Meteorological and Hydrological Institute
source: CMEMS BAL MFC NEMO model output converted to NetCDF
contact: servicedesk.cmems@mercator-ocean.eu
references: https://marine.copernicus.eu/
comment: Data on cropped native product grid. Horizontal velocities destagged
subset:source: ARCO data downloaded from the Marine Data Store using the MyOcean Data Portal
subset:productId: BALTICSEA_ANALYSISFORECAST_PHY_003_006
subset:datasetId: cmems_mod_bal_phy_anfc_P1M-m_202311
subset:date: 2024-06-11T06:25:27.439Z
st_tibble_new <-tidync(paste(covPath, "temperature_salinity","cmems_mod_bal_phy_anfc_P1M-m_1718087127439.nc",sep ="/")) |>hyper_tibble() |>mutate(date =as_datetime(time, origin ="1970-01-01")) |>mutate(month =month(date),day =day(date),year =year(date),month_year =paste(month, year, sep ="_") ) |>filter(year >2021)st_tibble <-bind_rows(st_tibble, st_tibble_new)# Loop through all year combos, extract the temperatures at the data locationsst_list <-list()for (i inunique(dat$month_year)) { d_sub <-filter(dat, month_year == i) st_tibble_sub <-filter(st_tibble, month_year == i)# Convert to rasterggplot(st_tibble_sub, aes(longitude, latitude)) +geom_point(size =0.1) st_raster <-as_spatraster(st_tibble_sub,xycols =3:4,crs ="WGS84", digits =3 )ggplot() +geom_spatraster(data = st_raster$bottomT, aes(fill = bottomT)) +scale_fill_viridis(option ="magma") +ggtitle(i)ggplot() +geom_spatraster(data = st_raster$sob, aes(fill = sob)) +scale_fill_viridis(option ="magma") +ggtitle(i)# Extract from raster d_sub$temp <- terra::extract( st_raster$bottomT, d_sub |> dplyr::select(lon, lat) )$bottomT d_sub$salinity <- terra::extract( st_raster$sob, d_sub |> dplyr::select(lon, lat) )$sob# Save st_list[[i]] <- d_sub}d_st <-bind_rows(st_list)